Magnet hypothesis: plant-pollinator interactions
Purpose: A test of the magnet hypothesis was examined in Mojave National Preserve by Ally Ruttan.
Hypothesis: Floral resource island created by shrubs and the associated beneficiary annual plants will positively and non-additively influence pollinator visitation rates.
Predictions:
(1) The frequency and duration of pollinator visitations to annuals is greater under shrubs than in the paired-open microsites (magnet H because of concentration).
(2) Annual plants under flowering entomophilous shrubs (Larrea tridentata) will have a higher frequency and duration of pollinator visitations than annual plants under anemophilous shrubs (Ambrosia dumosa) because of higher concentrations of suitable floral resources for pollinators (specificity of pollinator faciliation).
(3) Shrubs with annuals in their understory will have a higher frequency and duration of pollinator visitations than shrubs without annuals due to increased concentrations of floral resources for pollinators (reverse magnet effect and reciprocal benefits).
(4) Sites with both shrubs and annuals will have the highest frequency and duration of pollinator visitations to both the shrubs and the annuals (i.e. annuals under shrubs also with flowers are visited the most).
An interesting corollary is that there are appropriate floral resources for desert pollinators, that they discriminate, and that entomophilous and anemophilous shrubs facilitate flowering similarly.
Data wrangling
#libraries
library(tidyverse)
library(DT)
library(lubridate)
#meta-data
meta <- read_csv("data/meta-data.csv")
datatable(meta)#data
data.2015 <- read_csv("data/MNP.2015.csv")
data.2016 <- read_csv("data/MNP.2016.csv")
#merge
data <- rbind.data.frame(data.2015, data.2016)
data <- data %>% rename(net.treatment = treatment) %>% na.omit(data)
data$year <- as.character(data$year)
data$rep <- as.character(data$rep)
#tidy data to expand treatment column (current structure is a mix of three factors)
#microsite
data <- data %>% mutate(microsite = ifelse(net.treatment %in% c("SA", "SAA", "SX"), "Larrea", ifelse(net.treatment %in% c("OA"), "open", ifelse(net.treatment %in% c("AMB"), "Ambrosia", NA))))
length(unique(data$microsite))## [1] 3#annuals present
data <- data %>% mutate(annuals = ifelse(net.treatment %in% c("SA", "SAA"), "annuals", ifelse(net.treatment %in% c("OA"), "annuals", ifelse(net.treatment %in% c("AMB"), "annuals", "none"))))
length(unique(data$annuals))## [1] 2#target flowers
data <- data %>% mutate(target.flowers = ifelse(net.treatment %in% c("SA", "SX"), "shrub flowers", ifelse(net.treatment %in% c("AMB", "SAA", "OA"), "annual plant flowers", "NA")))
length(unique(data$annuals))## [1] 2#frequency counts in a separate dataframe (weighted by duration of recording)
frequency <- data %>% group_by(year, day, net.treatment, rep, microsite, annuals, target.flowers) %>% summarise(net.time = sum(total.duration), mean.visitation.duration = mean(visitation.duration), mean.floral.density = mean(floral.density), mean.insect.RTU = n_distinct(insect.RTU), count = n())
frequency$net.time <- as.numeric(frequency$net.time) #converts to total seconds
frequency$mean.visitation.duration <- as.numeric(frequency$mean.visitation.duration) #converts to total seconds
frequency$count <- as.numeric(frequency$count)
frequency$rate <- as.numeric(frequency$count)/frequency$net.time #weight by net time
frequency$proportion.visitations <- frequency$mean.visitation.duration/frequency$net.time
frequency$rate.per.flower <- frequency$count/frequency$mean.floral.density
datatable(frequency)#repeat above with RTU as factor
frequency.by.RTU <- data %>% group_by(year, day, net.treatment, insect.RTU, rep, microsite, annuals, target.flowers) %>% summarise(net.time = sum(total.duration), mean.visitation.duration = mean(visitation.duration), mean.floral.density = mean(floral.density), count = n())
frequency.by.RTU$net.time <- as.numeric(frequency.by.RTU$net.time) #converts to total seconds
frequency.by.RTU$mean.visitation.duration <- as.numeric(frequency.by.RTU$mean.visitation.duration) #converts to total seconds
frequency.by.RTU$count <- as.numeric(frequency.by.RTU$count)
frequency.by.RTU$rate <- as.numeric(frequency.by.RTU$count)/frequency.by.RTU$net.time #weight by net time
frequency.by.RTU$proportion.visitations <- frequency.by.RTU$mean.visitation.duration/frequency.by.RTU$net.time
frequency.by.RTU$rate.per.flower <- frequency.by.RTU$count/frequency.by.RTU$mean.floral.density
datatable(frequency.by.RTU)#repeat above but with RTU as species diversity response
#split out by year because of non-orthogonality
freq.2015 <- frequency %>% filter(year == 2015)
freq.2016 <- frequency %>% filter(year == 2016)Data visualization
#Ideal figures for publication.
#Two figures total: 1(a) rate, (b) duration and 2(a) rate per RTU and (b) mean visitation rate per RTU - maybe I think there are other options. ADD insect.RTU richness
#Higher-order treatment patterns in frequency####
#Collapsed single factor model
ggplot(frequency, aes(net.treatment, rate)) + geom_boxplot() + ylab("visitation rate") + facet_wrap(~year)
ggplot(frequency, aes(net.treatment, count)) + geom_boxplot() + ylab("visitations") + facet_wrap(~year)
ggplot(frequency, aes(net.treatment, rate.per.flower)) + geom_boxplot() + ylab("visitations per flower") + facet_wrap(~year)
ggplot(frequency, aes(net.treatment, proportion.visitations)) + geom_boxplot() + ylab("mean duration of visit per total duration recorded") + facet_wrap(~year)
ggplot(frequency, aes(net.treatment, mean.insect.RTU)) + geom_boxplot() + ylab("mean RTU richness") + facet_wrap(~year)
#net treatment per RTU
ggplot(frequency.by.RTU, aes(net.treatment, rate)) + geom_boxplot() + ylab("visitation rate") + facet_wrap(~insect.RTU*year)
ggplot(frequency.by.RTU, aes(net.treatment, proportion.visitations)) + geom_boxplot() + ylab("mean duration of visit per total duration recorded") + facet_wrap(~insect.RTU*year)
#Treatments separated
ggplot(frequency, aes(microsite, rate, fill = target.flowers)) + geom_boxplot() + facet_wrap(~year) + scale_fill_brewer(palette = "YlGn")
ggplot(frequency, aes(microsite, count, fill = target.flowers)) + geom_boxplot() + facet_wrap(~year) + scale_fill_brewer(palette = "YlGn")
ggplot(frequency, aes(microsite, rate.per.flower, fill = target.flowers)) + geom_boxplot() + facet_wrap(~year) + scale_fill_brewer(palette = "YlGn")
ggplot(frequency, aes(microsite, proportion.visitations, fill = target.flowers)) + geom_boxplot() + facet_wrap(~year) + scale_fill_brewer(palette = "YlGn")
ggplot(frequency, aes(microsite, mean.insect.RTU, fill = target.flowers)) + geom_boxplot() + facet_wrap(~year) + scale_fill_brewer(palette = "YlGn")
#relationships with sampling effort
ggplot(frequency, aes(net.time, count, color = year)) + geom_point()
ggplot(frequency, aes(net.time, count, color = year)) + geom_point() + facet_wrap(~microsite)
#floral density
ggplot(frequency, aes(mean.floral.density, rate, color = microsite)) + geom_point() + facet_wrap(~year)
ggplot(frequency, aes(mean.floral.density, proportion.visitations, color = microsite)) + geom_point() + facet_wrap(~year)
#relationships with sampling effort
ggplot(frequency, aes(net.time, rate, color = year)) + geom_point() + facet_wrap(~microsite*target.flowers)
ggplot(frequency, aes(net.time, mean.visitation.duration, color = year)) + geom_point() + facet_wrap(~microsite*target.flowers)
EDA
#test distributions and explore outliers
summary(frequency)## year day net.treatment
## Length:266 Length:266 Length:266
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## rep microsite annuals
## Length:266 Length:266 Length:266
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## target.flowers net.time mean.visitation.duration
## Length:266 Min. : 300 Min. : 0.000
## Class :character 1st Qu.: 4500 1st Qu.: 1.475
## Mode :character Median : 23010 Median : 26.248
## Mean : 64855 Mean : 58.067
## 3rd Qu.:100660 3rd Qu.: 47.584
## Max. :542929 Max. :2484.114
## mean.floral.density mean.insect.RTU count rate
## Min. : 5.00 Min. :1.000 Min. : 1.00 Min. :0.0001774
## 1st Qu.: 23.13 1st Qu.:2.000 1st Qu.: 3.00 1st Qu.:0.0002130
## Median : 71.26 Median :3.000 Median : 8.00 Median :0.0002780
## Mean : 98.74 Mean :3.342 Mean : 15.14 Mean :0.0005540
## 3rd Qu.:200.00 3rd Qu.:4.000 3rd Qu.: 22.00 3rd Qu.:0.0011111
## Max. :200.00 Max. :9.000 Max. :109.00 Max. :0.0033333
## proportion.visitations rate.per.flower
## Min. :0.000e+00 Min. :0.00500
## 1st Qu.:9.565e-05 1st Qu.:0.02134
## Median :4.459e-04 Median :0.10158
## Mean :3.221e-03 Mean :0.62738
## 3rd Qu.:1.633e-03 3rd Qu.:0.86012
## Max. :2.117e-01 Max. :6.05000require(fitdistrplus)
descdist(freq.2015$rate, boot = 1000)
## summary statistics
## ------
## min: 0.0001873882 max: 0.001337793
## median: 0.0003267998
## mean: 0.0004887686
## estimated sd: 0.0003435598
## estimated skewness: 1.110484
## estimated kurtosis: 2.810992descdist(freq.2016$rate, boot = 1000)
## summary statistics
## ------
## min: 0.0001773679 max: 0.003333333
## median: 0.0002480947
## mean: 0.0006144842
## estimated sd: 0.000504536
## estimated skewness: 1.331961
## estimated kurtosis: 7.146892detach("package:fitdistrplus", unload = TRUE)
#explore temperature on count and mean visitation ratesModels
#GLM for count and weight by net.time (alt - use MASS and glm.nb)
#library(MASS) #need for glm.nb
#all codes
#2015 counts
m <- glm(count~net.treatment + mean.floral.density, family = "poisson", weight = net.time, data = freq.2015)
anova(m, test = "Chisq") ## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 127 47238655
## net.treatment 3 13736996 124 33501659 < 2.2e-16 ***
## mean.floral.density 1 13630462 123 19871197 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")## $lsmeans
## net.treatment lsmean SE df asymp.LCL asymp.UCL
## OA 2.045110 0.0004169528 NA 2.044293 2.045927
## SA 2.614917 0.0005830291 NA 2.613774 2.616060
## SAA 1.736211 0.0004255844 NA 1.735377 1.737045
## SX 1.354487 0.0007356275 NA 1.353045 1.355929
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## OA - SA -0.5698067 0.0007644926 NA -745.340 <.0001
## OA - SAA 0.3088991 0.0002022540 NA 1527.283 <.0001
## OA - SX 0.6906229 0.0007383659 NA 935.340 <.0001
## SA - SAA 0.8787058 0.0007683106 NA 1143.686 <.0001
## SA - SX 1.2604296 0.0009583930 NA 1315.149 <.0001
## SAA - SX 0.3817238 0.0007455663 NA 511.992 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimates#2016 counts
m <- glm(count~net.treatment + mean.floral.density, family = "poisson", weight = net.time, data = freq.2016)
anova(m, test = "Chisq") ## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 137 174945050
## net.treatment 4 15411680 133 159533369 < 2.2e-16 ***
## mean.floral.density 1 8872062 132 150661307 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")## $lsmeans
## net.treatment lsmean SE df asymp.LCL asymp.UCL
## AMB 5.2074833 0.0005064262 NA 5.2064907 5.2084759
## OA 5.2481438 0.0004770340 NA 5.2472089 5.2490788
## SA -1.0445768 0.0023971369 NA -1.0492751 -1.0398785
## SAA 5.3522565 0.0004641314 NA 5.3513468 5.3531661
## SX -0.7880837 0.0019263237 NA -0.7918592 -0.7843082
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## AMB - OA -0.04066053 1.176853e-04 NA -345.502 <.0001
## AMB - SA 6.25206012 2.597272e-03 NA 2407.164 <.0001
## AMB - SAA -0.14477316 1.130315e-04 NA -1280.821 <.0001
## AMB - SX 5.99556701 2.170318e-03 NA 2762.529 <.0001
## OA - SA 6.29272065 2.583621e-03 NA 2435.621 <.0001
## OA - SAA -0.10411263 9.853844e-05 NA -1056.569 <.0001
## OA - SX 6.03622754 2.153962e-03 NA 2802.383 <.0001
## SA - SAA -6.39683328 2.578069e-03 NA -2481.250 <.0001
## SA - SX -0.25649310 2.889407e-03 NA -88.770 <.0001
## SAA - SX 6.14034017 2.147300e-03 NA 2859.563 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 5 estimates#repeat above for mean.visitation.duration
#2015
m <- glm(mean.visitation.duration~net.treatment + mean.floral.density, family = "gaussian", weight = net.time, data = freq.2015)
anova(m, test = "Chisq") ## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: mean.visitation.duration
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 127 7.3579e+10
## net.treatment 3 1.1426e+10 124 6.2153e+10 4.864e-05 ***
## mean.floral.density 1 1.5637e+06 123 6.2151e+10 0.9556
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")## $lsmeans
## net.treatment lsmean SE df asymp.LCL asymp.UCL
## OA 71.412514 29.25204 NA 14.07957 128.74546
## SA 5.832637 40.04779 NA -72.65959 84.32486
## SAA 145.904619 29.89835 NA 87.30492 204.50431
## SX 12.115824 39.10838 NA -64.53519 88.76684
##
## Results are given on the identity (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## OA - SA 65.579877 56.58232 NA 1.159 0.6528
## OA - SAA -74.492105 21.63099 NA -3.444 0.0032
## OA - SX 59.296690 47.54971 NA 1.247 0.5968
## SA - SAA -140.071982 56.69949 NA -2.470 0.0646
## SA - SX -6.283187 56.59297 NA -0.111 0.9995
## SAA - SX 133.788795 47.99355 NA 2.788 0.0272
##
## P value adjustment: tukey method for comparing a family of 4 estimates#2016
m <- glm(mean.visitation.duration~net.treatment + mean.floral.density, family = "gaussian", weight = net.time, data = freq.2016)
anova(m, test = "Chisq") ## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: mean.visitation.duration
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 137 1.3624e+12
## net.treatment 4 1.0011e+11 133 1.2623e+12 0.03313 *
## mean.floral.density 1 8.5266e+08 132 1.2614e+12 0.76516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#posthoc test
require(lsmeans)
lsmeans(m, pairwise~net.treatment, adjust="tukey")## $lsmeans
## net.treatment lsmean SE df asymp.LCL asymp.UCL
## AMB -60.98750 348.6464 NA -744.3219 622.3469
## OA 136.40291 330.4878 NA -511.3413 784.1471
## SA 167.37093 649.1859 NA -1105.0100 1439.7519
## SAA -59.88931 322.6658 NA -692.3026 572.5240
## SX 183.58414 622.0679 NA -1035.6465 1402.8148
##
## Results are given on the identity (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## AMB - OA -197.390411 75.65485 NA -2.609 0.0687
## AMB - SA -228.358432 945.87912 NA -0.241 0.9993
## AMB - SAA -1.098192 73.96065 NA -0.015 1.0000
## AMB - SX -244.571636 927.47694 NA -0.264 0.9989
## OA - SA -30.968021 929.59326 NA -0.033 1.0000
## OA - SAA 196.292219 67.38350 NA 2.913 0.0295
## OA - SX -47.181225 910.86211 NA -0.052 1.0000
## SA - SAA 227.260240 922.95167 NA 0.246 0.9992
## SA - SX -16.213203 535.98700 NA -0.030 1.0000
## SAA - SX -243.473443 904.08293 NA -0.269 0.9988
##
## P value adjustment: tukey method for comparing a family of 5 estimates#treatments split out
#how to test?Interpretation
Net.treatment models
1. Using the simplest models with all treatment levels aggregated but split by years, larrea is the most important magnet for frequency of visits (weighting for net.time recorded and using floral density as a covariate). This is a reasonable test because it captures enough of the variation and address non-orthogonality levels.
2. Mean visitation rate in 2015 significantly differed between larea with flowers and open and larea with flowers and annuals versus larrea without. This strongly suggests a double-magnet effect this season. In 2016, larrea flowering with annuals was also significantly greater than open microsites with annuals.
Double-magnet supported for frequency too.
Then check each prediction. Need to think over how to handle double-magnet H.